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ABSTRACT 

The effects of dust extinction on the departure from axi- symmetry in the 
morphology of planetary nebulae (PNs) are investigated through a comparison 
of the radio free-free emission and hydrogen recombination line images. The dust 
exctinction maps of five compact PNs are derived using high-resolution (~0'.'l) 
Ha and radio images maps of the HST and VLA. These extinction maps are 
then analyzed by an ellipsoidal shell ionization model including the effects of 
dust extinction to infer the nebulae's intrinsic structure and orientation in the 
sky. This method provides a quantitative analysis of the morphological structure 
of PNs and represents a step beyond qualitative classification of morphological 
types of PNs. 

Subject headings: ISM: planetary nebulae: general — stars: AGB and post-AGB 
— stars: supernovae: general 
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1. Introduction 

In general, PNs have axisymmetric structures. The classification of PNs into three basic 
types, round, elliptical and bipolar (Balick 1987), is based on their axisymmetric appear- 
ances. However, many PNs also possess non-axisymmetric structures. Soker & Hadar (2002) 
proposed a scheme to classify PNs according to their departure from axisymmetry. They 
considered six categories: (1) PNs where the central star is not at the center of the nebula; 
(2) PNs having one side brighter than the other; (3) PNs having unequal size or shape of 
the two sides; (4) PNs where the symmetry axis is bent, e.g., the two lobes in a bipolar 
PN are bent toward the same side; (5) PNs where the main departure from axisymmetry 
is in the outer regions, e.g., an outer arc; (6) PNs showing no large-scale departure from 
axisymmetry. They found that about 50% of the analyzed sample has large-scale departures 
from axisymmetry. 

The problem of departure from axisymmetry in PNs is usually investigated through the 
examination of optical images. This poses a problem because dust also plays a role in making 
such asymmetry. For example, in the sample analyzed by Soker & Hadar (2002), more than 
half of the PNs they found to have departure from axisymmetry are the types with unequal 
intensity or size. This kind of asymmetric appearance can also be caused by dust extinction, 
even if the intrinsic brightness of the nebula is axisymmetric. It is very difficult to draw 
conclusions about the real shapes of PNs based only on their apparent optical forms. In 
order to remove the bias caused by dust extinction while investigating the shapes of PNs, it 
is necessary to understand the dust distribution in PNs. 

Asymptotic Giant Branch (AGB) stars eject large amount of dust and gas in the form 
of stellar winds. These stellar grains is a major source of dust grains in the Galaxy (Gehrz 
1989). Although the details of dust formation are still poorly understood, the properties 
of circumstellar dust have been extensively studied in the infrared, e.g., by the Infrared 
Astronomy Satellite (IRAS) and Infrared Space Observatory (ISO) missions. Since PNs 
descend from AGB stars, the remnants of the AGB dust envelopes should still be present 
in PNs (Kwok 1982). The presence of dust in PNs was first detected in NGC 7027 via its 
infrared emission (Gillett, Low & Stein 1967). This infrared excess peaks around 30 /iin and 
has a color temperature of ~100 K. Although the existence of dust in PNs first came as a 
surprise, its widespread presence is now confirmed by IRAS observations (Pottasch et al. 
1984; Zhang & Kwok 1991). 

While the presence of dust in PNs is well established by infrared spectroscopy, very 
little is known about its spatial distribution. Recent development in mid-infrared cameras 
has just begun to explore the dust distribution of PNs (see, e.g., Volk & Kwok 2003; Su 
et al. 2004). However, a full mapping of the dust distribution of PNs has to wait until the 
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development of large-format mid-infrared cameras with high resolution and high sensitivity. 
An alternate and underutilized method to derive the dust distribution is by comparing the 
hydrogen recombination line map of a PN with its corresponding radio free-free continuum 
map. The recombination line (e.g. Ha) and free-free continuum emission are independent 
tracers of the distribution of the ionized gas, and both of their fluxes are proportional to 
n 2 e V under optically thin conditions, where n e is the electron density and Vis the volume of 
the nebula. If dust is present, they will have no effect on the long wavelength radio emission, 
but the Ha emission will be attenuated by dust. Therefore the dust distribution can be 
indirectly obtained by comparing the Ha and radio maps, assuming the dust is associated 
with the ionized gas. This method has been applied to different objects by various authors 
(e.g., Vazquez et al. 1999; Eyres et al. 2001). In this paper, we have derived optical depth 
maps of five compact PNs obtained with the same technique. Compact PNs are selected 
for this study because their small angular sizes allow the assumption that any foreground 
interstellar extinction is homogeneous over the entire object, so that the different degrees of 
extinction over a PN are only related to its internal dust distribution. 



2. Derivation of the optical depth distribution 

In a fully ionized region, both the free-free emission and the recombination fluxes are 
proportional to n%V, where n e is the electron density. The ratio of the free-free flux to the 
H/3 recombination line fluxes under Case B (where recombination to the ground state is 
excluded) is given by the following equation (Pottasch 1984): 

^ (ff) - 2.51 x lO^u'^Y Jan8ky 



F(K/3) ' e ergcm^s- 1 

where T e is the electron temperature, and 

n(He+) + 3 7 n(He++) 



Y = l + -^ ^ + 3.7— - (2) 

Tip Tip 

corrects for the contribution to the free-free flux from He ions. For a He to H number ratio 
of 0.11 and for He to be equally divided between singly and doubly ionized forms, Y has a 
value of 1.258. Under Case B, the effective recombination coefficients for Ha and H/3 have a 
ratio of 2.85 (Hummer & Storey 1987). For T e = 10 4 K, eq. 1 can be written as 

F exp (Ha) = 6.85 x 10" 10 (z//GHz) ai ( 1 S iy /Jy) erg cm" 2 s" 1 , (3) 

which can be used to calculate the expected flux of Ha emission without dust extinction 
from the flux density of radio continuum emission. 
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In order to derive the extinction map, we begin with the radio continuum and Ha maps 
of the same PN at similar angular resolutions. The expected Ha emission at each pixel is 
derived from the radio continuum map, and then divided by the observed Ha emission map 
for each object to obtain an extinction map, in which the dust optical depth at each pixel 
can be calculated. The intensity emitted by an object and received by the observer is related 
by the optical depth r as 



obs 



-^expC ? (4) 



where F exp is the intensity emitted by the object, i.e., the intensity expected to be received 
by the observer if no dust is present. Therefore the dust optical depth for Ha emission is 
obtained by taking the natural logarithm of the ratio of the expected and observed Ha flux: 



3. Sample selection and data processing 

In this paper, we have selected five compact PNs with known high infrared excesses. 
Figures 1 show the spectral energy distributions (SEDs) of these five PNs from optical to 
radio wavelengths, illustrating the IR excess caused by dust. The ground-based optical 
data are from Acker et al. (1992) and Acker, Marcout & Ochsenbein (1981). The near- 
infrared J, H, and if-band photometry are from the Two Micron All Sky Survey (2MASS), 
extracted from the 2MASS All-Sky Catalog of Point Sources (Cutri et al. 2003). The mid- 
infrared photometry includes IRAS photometry, IRAS low resolution spectra (LRS), and 
recent measurements from the Midcourse Space Experiment (MSX). The MSX measurements 
were extracted from the MSX Point Source Catalog (Egan et al. 2003). Radio photometric 
measurements (mostly at A 2 and 6 cm) are from Aaquist & Kwok (1990), Basart & Daub 
(1987), and this study. As clearly shown in Fig. 1, the infrared excess of the PNs is well 
represented by a single blackbody curve with a temperature ~ 200 K. 

High resolution Ha and radio images of these PNs are available from the HST and the 
VLA respectively. The basic information for the five objects is listed in Table 1. With 
these high resolution images, the point-to-point spatial extinction distribution in PNs can 
be obtained with angular resolutions as high as ~ 0.1 arcsec. 

For IC 5117, M 1-61, M 2-43, and M 3-35, the HST WFPC2 Ha images were obtained 
under the GO program 8307 (Kwok, Su & Sahai 2003). For BD+30°3639, the HST WFPC2 
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Ha image was first presented in Harrington et al. (1997) and the calibrated images were 
obtained from the HST archive. The only image processing needed was to remove cosmic 
rays with the task CRREJ in IRAF. The angular resolution FWHM of these images ranges 
from 0"085 to O'/ll. 

For BD+30°3639, the A6 cm radio image was from Bryce et al. (1997). It was obtained 
by combining VLA and Multi-Element Radio Linked Interferometer Network (MERLIN) 
observations, resulting in a very high angular resolution of ~ 0'/08. For the remaining four 
objects, the radio continuum maps were observed by Aaquist & Kwok with the VLA at A2 
cm in A-configuration in March 1998. The unpublished data were retrieved from the VLA 
archive and reduced using the package AIPS following the standard procedures of continuum 
data reduction, as summarized in Appendix A of the AIPS Cookbook (NRAO 2003). After 
the data were calibrated, the final CLEAN images were generated with robust weighting 
in order to obtain images with minimum-noise levels while maintaining the highest angular 
resolution. This ranges from ~ 0'/12 to ~ 0'.'15 in the final CLEAN maps. 

The operations taken to combine the radio and optical data and to derive the optical 
depth map are performed in the following sequence. The HST images were first oriented 
north-up with Karma software. Subsequent regridding and other steps were performed in 
AIPS. Because the coordinates of the i^ST images have a uncertainty of ~ 0'/5, there is usually 
an offset in the positions of the same object in the radio and optical images. To correct this, 
the Ha images were shifted and aligned by eye to match with the radio continuum maps 
using features common in both images. Table 3 gives the offset of the coordinates for each 
object; a positive sign indicates the optical image has been shifted eastward or northward. 
Each object's Ha image was convolved with an extended Gaussian to match the beam size 
of the corresponding radio image, except for BD+30°3639. Because the radio continuum 
image of BD+30°3639 has a higher resolution than its Ha image, it was convolved with an 
extended Gaussian to match the FWHM of the Ha image. The radio image was transformed 
so that its pixel coordinate grid was consistent with that of the HST image. Division and 
natural logarithm operations were then performed to obtain an optical depth map as in eq. 5. 



4. Results and Analysis 

Extinction maps showing the dust optical depth for the five compact PNs are presented 
in Fig. 2. For comparison, the Ha brightness and extinction maps for each object are shown 
side by side, and the corresponding radio continuum contours are superposed on both maps. 
All images are oriented with north up and east to the left. 
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4.1. BD+30°3639 



The radio continuum image of BD+30°3639 shows a ring-like structure with two peaks 
of roughly the same brightness in the north and south regions. The Ha image also shows a 
ring-like structure with two peaks. However, the ring shows unequal brightness in the east 
and west regions, which indicates the ring suffers more dust extinction in the west as shown 
in the optical depth map. The two high-extinction regions in the northeast and southwest 
are in roughly the same positions as the pair of high-velocity molecular knots found with 
CO line imaging by Bachiller et al. (2000). 



The radio continuum image of IC 5117 shows a shell structure with two emission peaks. 
Such a shape can be simply explained as a prolate ellipsoidal shell projected onto the plane 
of the sky (Aaquist & Kwok 1991). The Ha image shows a similar structure with two pairs 
of faint bipolar lobes oriented at different angles (see image displayed by Kwok, Su & Sahai 
2003). The lobes are not seen in the radio continuum because of relatively low sensitivity 
for the radio map. The extinction map shows a similar structure to that of the radio map, 
where the high-extinction region follows the region of high brightness in radio. However, the 
two extinction peaks fall slightly outside of the radio peaks. A similar situation has been 
found in the dust extinction map of NGC 7027 (Walton et al. 1988). 



The radio continuum image shows a shell structure with two pronounced emission peaks 
of roughly the same brightness. In the Ha image, the eastern peak is weaker than the western 
peak, and there is at least one pair of faint bipolar lobes in the outer region (see image 
displayed by Kwok, Su & Sahai 2003). The extinction map shows the optical depth is higher 
in the region of brighter radio emission, with a strong extinction peak near the eastern radio 
peak. 



4.2. 



IC 5117 



4.3. 



M 1-61 



4.4. 



M 2-43 



The radio continuum image shows a shell structure, with the southern peak brighter 
than the northern one. The Ha image shows similar structure with slightly extended faint 
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emission outside the radio emission region. However, the northern peak of Ha emission is 
brighter than the southern one. As a result, although the structure of the extinction map 
mainly follows that of the radio map, the optical depth in the southern part is significant 
higher than the northern part. The extinction peaks also falls slightly outside of the radio 
peaks. 



The radio continuum image shows a deep central emission minimum and two surround- 
ing peaks with concave-like contours at roughly the same brightness. The Ha image shows 
two peaks at approximately the same location with two pairs of faint bipolar lobes outside 
(see image displayed by Kwok, Su & Sahai 2003). The extinction map has similar structure 
to that of the radio map, with the high-extinction region coinciding with the concave-like 
region. 



Figure 3 show the uncertainties propagated from the radio and optical images. The 
uncertainty associated with the optical depth r is given by 



where o" ra dio and (Tho are the rms noise of the radio continuum and optical Ha images, 
respectively. Since the optical images have higher signal-to-noise ratios, the uncertainty of 
the optical depth is dominated by the radio image. Hence the uncertainty map of the optical 
depth and the radio continuum map have similar structure. The optical depth uncertainty 
ranges from At ~ 0.01 to 0.4, where smaller values lie in the regions with stronger radio 
emission. 

The calculations of optical depth above use several assumptions. To investigate whether 
the optical depth is sensitive to these assumptions, the uncertainties propagated from them 
should also be examined. Among these approximations, the electron temperature T e gives 
the most significant uncertainty. While most PNs have T e in the vicinity of 10,000 K, the 
electron temperatures of PNs cover a range of 5,000 - 15,000 K, with T e ~ 20,000 K in some 
extreme cases. As a result, different choices of T e will give different values of optical depth. 



4.5. 



M 3-35 



4.6. Uncertainties in the Optical Depth 





To calculate the uncertainty propagated from the assumption that T e = 10, 000 K, we 
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estimate the uncertainty for T e to be roughly AT e = ±5,000 K. Since r oc lnF cxp (Ha) oc 
lnT" 53 , the uncertainty propagated into the optical depth is At oc 0.53 x Ap^ oc 0.53 x 
io'ooo ~ 0.27. This value is an order of magnitude smaller than the derived optical depth. 
This uncertainty from the electron temperature will only shift the values of the derived 
optical depths by a constant, but will not change the overall appearance of the extinction 
caused by its dust in a PN. 



4.7. Dust distribution 

All five PNs in this study show that the Ha emission peaks and the radio continuum 
emission peaks are in the same regions. In four of them, the peaks of the dust extinction 
also reside in the same regions. These results suggest that the ionized gas and the dust have 
similar distributions, so at least part of the dust is mixed with the ionized gas. 

On the quantative scale, the faintest parts of the optical image correlate with the peaks 
of the extinction distribution, as is expected. As an independent check, Kastner et al. (2002) 
compared the X-ray emission map of BD+30°3639 with its extinction map derived from 
optical and infrared images, and showed that the peak of X-ray emission also lies in the 
region with minimum extinction. 

We note that the radio images, Ha images, and the dust optical depth maps only rep- 
resent the column density distribution along the line of sight. In order to remove projection 
effects and reconstruct the three dimensional ionized gas and dust distributions in PNs, we 
need to develop a 3-D model. The details of the model analysis are presented in the next 
section. 



5. The Three Dimensional Dusty Ellipsoidal Shell (DES) Model 

While the classification of PN morphology is based on appearance, this approach is vul- 
nerable to projection effects. For example, elliptical and bipolar nebulae can appear round 
when the nebulae are viewed pole-on. In order to study their intrinsic properties, there have 
been many attempts to find the detailed three-dimensional structures of PNs (Khromov & 
Kohoutek 1968). Since many PNs have the appearance of limb- brightened elliptical rings 
that are brightest on their minor axes, it is natural to consider models in which the nebulae 
have an elongated, hollow shell structure. Models with this simple geometry, such as the 
prolate ellipsoidal shell discussed by Masson (1989, 1990), are successful in reproducing many 
of the shapes observed in PNs. This was followed by the work of Aaquist & Kwok (1996) 
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who developed a prolate ellipsoidal shell (PES) model to fit the radio images of PNs. Their 
model contains a spherical shell with both radial and latitude-dependent density gradients, 
and the shell is ionized by a central star to various depths at different latitudual directions. 
The ionized shell, when viewed at different angles, results in a variety of morphologies seen 
in PNs. The parameters derived from the PES model reveal the asymmetric distribution of 
the ionized gas in the nebula. The PES model was later called the ES model and explored 
extensively by Zhang & Kwok (1998) to investigate the statistical distributions of the asym- 
metry parameters. However, because of a typographical error in the source code, there are 
some uncertainties in the derived model parameters in Zhang & Kwok (1998). 

In order to simulate dust extinction effects, we have modified the ES model to include a 
dust distribution with a constant gas-to-dust ratio, and we refer to this model as the dusty 
ellipsoidal shell (DES) model. The DES model has been used to analyze the observational 
images, including radio continuum, Ha, and dust extinction maps to constrain the ionized 
gas and dust distributions in compact PNs. These different images of the same object provide 
a way to study the distribution of dust in PNs, and also give better constraints on the model 
parameters than a single radio or optical image would. One advantage of using the model is 
that it helps to remove projection effects and to construct the three-dimensional structures 
of PNs. The DES model will give a quantitative classification of PNs morphology rather 
than the standard descriptive terms such as round, elliptical, or bipolar. 



5.1. The Geometry of the DES model 

The geometry of the DES model is shown in Figure 4. A single (hot) star is located inside 
a shell of ionized gas. The inner region of the nebula is assumed to be empty, supposedly 
swept clean by a wind from the central star. For simplicity, the empty region is assumed 
to be a prolate or oblate spheroid, depending on the choice of polar radius a and equatorial 
radius b. The shape of the shell is the result of the interacting winds process, where either 
the slow AGB wind is equatorially enhanced, or the fast central star wind is non-isotropic 
(e.g. collimated), or both. The distance from the star to the inner edge of the gas R s can be 
calculated from a, b and the polar angle 9. The UV photons from the central star ionize the 
surrounding gas out to the distance Ri, which depends upon R s and the density structure 
of the surrounding nebula. 

The model assumes the central star is emitting ionizing photons isotropically, and the 
ionized nebula is ionization-bounded along any radial directions from the central star. The 
depth of the ionization front along any radial vector is given approximately by 
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dL f Ri 



n 2 (r, 9, 4>)r 2 dr , 



(7) 



where dL/dQ, the number of local ionizing photons (L) per unit solid angle (Q), is set equal 
to the total number of recombinations within the local ionized region. In this equation, «g is 
the recombination coefficient excluding captures to the ground level (Case B approximation), 
n(r, 9, 0) is the density structure of the swept-up shell, and r, 9, are the spherical coordi- 
nates. Ri is a function of 9 and that depends on the density distribution. For simplicity, a 
separable density law is assumed, so that n(r, 9, 0) = norj r rjeV<t>] here no is the density at the 
inner edge, where (r, 9, 0) = (R s , 0, 0), and i] r , i] e , r]^ are dimensionless quantities describing 
the density relative to n : 



ir/Rs)- 



(8) 



Ve = fe x < 



V<p = U ■ 
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2tt-2# x ' 



if < < 



7T 



7T 



7T 



if - < 9 < 7T 

2 - _ 



(9) 



(10) 



In this formulation, 7 gives the logarithmic radial density gradient, (5 gives the pole-to 
equator density ratio, a governs the shape of the angular density function, and fe,f<j, are 
additional asymmetry factors: 



fe = 5 e + (1 - 5 e )exp[cos(9 - 9 ) - 1}; < 5 e < 1 
U = $4> + (! _ ^)exp[cos(0 - </> ) - 1]; 0<^<1, 



(11) 
(12) 



where these control asymmetry in the latitude (#) and longitude (0) directions, respectively. 
These two asymmetry factors are used to produce the departure from axisymmetry seen in 
many PNs. For simplicity, the phases 9 and O are both assumed to be 0°, so that fe, f<p 
are maximized (=1) for = 0° and = 0°, and minimized for = 180° and = 180°. 

With the given radial density law r/ r , eq. 7 can be integrated to give a simple analytical 
solution: 

^ = olb f {n r] r r]e'n < j > ) 2 r 2 dr 



dn 



Ri 



(n r}eV<t>) 2a B / r 2l+2 R 2l dr 
Rs 
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= (noVeVt) ®bR s _^ - 3 • (13) 

Once the parameters describing the inner boundary, the density law (7, rje, rj^), and the 
equatorial axis shell thickness Ab = [Ri — R s ] g=90 o are chosen, the outer ionized boundary 
Ri(9, 4>) can be calculated from eq. 13. 

The ionized shell is then projected onto the plane of the sky at different inclination 
angles i. For % = 0°, the nebula is seen pole-on, i.e., the line of sight is perpendicular to 
the equatorial plane. The edge-on case has % = 90°, with the line of sight in the equatorial 
plane. The projected emission from the ionized region is calculated by assuming that the 
intensity of the recombination Ha line and the radio continuum emission is proportional to 
the integration of the square of the density function along the line of sight (eqs. 8 - 12). The 
calculations were done inside a cubic box containing the nebula. 

To simulate optical extinction effect, we assume that the dust is well mixed with the gas 
in the ionized region. Thus, a constant gas-to-dust ratio is assigned throughout the ionized 
region, i.e., the amount of extinction is directly proportional to the density function at each 
integration point. The Ha emission is then attenuated by the dust extinction using eq. 4. 



5.2. The Model Images 

The DES model allows us to explore how the morphology of a nebula changes with 
various parameters, especially a and (3. A series of 128 x 128 pixel model images has been 
produced with fixed parameters of a = 20, b = 18, and Ab = 6. The model images presented 
here have no radial density gradient, i.e., 7 = 0. This is a fairly good approximation, because 
most swept-up shells are thought to be relatively thin, allowing little density change with 
radius. We have computed the projected images for inclination angles % = 0°,30°,60°, and 
90° with different combinations of (3 and a as shown in Figures 5-8. 

Figure 5 shows the model images with emission calculated from the integration of the 
square of the density function along the line of sight. This simulates the radio continuum 
images of the long- wavelength free- free emission not affected by dust. 

Figure 6 displays the input dust distribution. The images are computed from the inte- 
gration of the amount of extinction along the line of sight all the way through the nebula. 
Because the gas-to-dust ratio is constant inside the ionized region, the dust distributions 
appear similar to the model radio images in Figure 5. If the dust temperature is uniform, 
and the cloud is optically thin for the dust thermal emission, the flux emitted by the dust 
is directly proportional to the absorption coefficient. Therefore, these dust distributions can 
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mimic infrared images of the dust emission. 

Figure 7 gives the model images when the dust extinction is taken into account. As in 
Figure 5, emission is calculated from the integration of the square of the electron density, 
but now it is attenuated by the dust in front of it along the line of sight. This simulates 
the optical Ha images, which, unlike the radio continuum maps, are affected by dust. Some 
images clearly show an asymmetric appearance caused by dust extinction (e.g., the image of 
a = 1, (3 = 0.1, and % = 30°). 

Figure 8 displays output extinction optical depth maps, calculated as the natural loga- 
rithm of the ratio of the model radio and optical images (eq. 5). Because only dust in front 
of the emission can attenuate the light, the amount of dust in the extinction maps (Fig. 8) 
is less than the actual amount of dust in the nebulae (Fig. 6). Comparing these images with 
the input dust distributions in Figure 6, the amount of missing dust for different inclination 
angles and morphologies can be estimated. By comparing the radio and Ha images, the neb- 
ular geometry can be estimated and the asymmetry effect caused by dust for this viewing 
geometry can be removed and reveal any intrinsic asymmetry of the gas distribution. 

6. Comparison of Observed and Model Images 

In order to constrain the parameters that describe the geometry and angular density 
distribution of the PNs presented, a set of model images were made to compare to each indi- 
vidual nebula's observed radio continuum and Ha images, extinction maps, and in one case 
(BD+30°3639), an infrared dust emission image (from Volk & Kwok 2003). The comparison 
was made by visual inspection. The faint emission from the outer parts of each nebula was 
chosen to match the observed Ha images because their S/N ratios are much higher than 
those for the radio continuum images. For each object, a set of parameters was found that 
generate model images matching most of the observed images. The images comparing the 
best models to the observations for these five nebulae are shown in Figures 9 - 13. The 
parameters used to generate the model images are listed in Table 8. 

6.1. BD+30°3639 

The model of BD+30°3639 has a relatively small equatorial thickness Ab compared to 
the size of the inner empty region, which has a prolate shape with an axial ratio of a/b = 1.2. 
The small value of (3(= 0.1) suggests that it has a large equatorial to polar density ratio. 
A small asymmetry 5<f, = 0.9 in the longitude direction is added to explain the unequal 
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intensity of the two sides observed in the radio image. Dust well-mixed with ionized gas is 
able to reproduce the anti-correlation of the optical image and extinction map. Due to its 
low inclination angle (i = 15°), it has a round appearance. This low inclination is consistent 
with the value of % ~ 10° found by Li, Harrington & Borkowski (2002), who also fit the radio 
image and the spectra by approximating the nebula as an ellipsoid. But a comparison with 
the model image in Fig. 7 (second column from the right) suggests that BD+30°3639 would 
appear as a prominent bipolar nebula if viewed edge-on. 

6.2. IC 5117 

The model of IC 5117 also has a relatively small equatorial thickness compared to the 
inner empty region, with a prolate axial ratio a/b = 1.67. It has a small pole-to-equator 
ratio of /3 — 0.1. No asymmetry factor in the longitude direction is needed. The brighter 
region toward the east in the optical image and the higher extinction toward the west can 
be reproduced by a dust distribution with a constant gas-to-dust ratio. The model images 
have an intermediate inclination angle of i — 40°. 

6.3. M 1-61 

The model of M 1-61 has an equatorial thickness comparable to the dimension of the 
inner empty region, which has a prolate axial ratio of a/b = 1.2. It has a small pole-to- 
equator ratio of (5 — 0.2. A small asymmetry 5$ = 0.9 in the longitude direction is needed 
to explain the mildly brighter western peak in the radio image. The anti-correlation of the 
optical bright peak in the west and the strong extinction peak in the east can be reproduced 
by the well-mixed dust. However, it is not as prominent as the observed images. The model 
images has a high inclination angle of % = 70°, and they suggest that a higher dynamic range 
image of M 1-61 will reveal a pair of bipolar lobes. 

6.4. M 2-43 

The model of M 2-43 also has a relatively small equatorial thickness compared to the 
inner empty region with a prolate axial ratio a/b = 1.67. It has a relatively small pole-to- 
equator ratio of f3 — 0.3. A mild asymmetry 5$ = 0.7 in the longitude direction is added 
to reproduce the brighter southern peak in the radio image. Although a simple constant 
gas-to-dust ratio is able to generate an anti-correlation between the Ha emission and optical 
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depth, it cannot reproduce optical and extinction maps that match the observed images well. 
The model images have an intermediate inclination angle of % — 45°. 

6.5. M 3-35 

The model of M 3-35 has an equatorial thickness comparable to the dimension of the 
inner empty region, which has an oblate axial ratio of a/b = 0.83. This is the only object 
that requires an oblate model. It has a intermediate pole-to-equator ratio of f3 — 0.4. No 
asymmetry factor in the longitude direction is needed to reproduce the radio image. No 
obvious anti-correlation between optical and extinction maps can be found. The model 
images are highly inclined to the line of sight, with i = 80°. As in M 1-61, a higher dynamic 
range image is expected to reveal the bipolar lobes in the NE-SW directions. 

7. Discussion 

As we can see from Figures 9 - 13, the model images can generally simulate the observed 
images. The remaining differences could be accounted for by several effects. For example, 
the assumed model geometry may be too simple to reproduce all of the complex small-scale 
structure of PNs. Other causes could include clumpy dust distributions that do not follow 
a simple constant gas-to-dust ratio. Three out of five PNs show an asymmetric density dis- 
tribution in the longitudinal direction, a departure from axisymmetry. As summarized by 
Soker & Rappaport (2001), four main processes can result in the departure from axisymme- 
try in PNs: interaction with the interstellar medium, mass-loss enhancements due to cool 
starspots, or a wide or a close binary companion. They investigated the effects of companions 
and estimated that about 27% of all PNs are expected to acquire non-axisymmetric structure 
from binary interactions. Nevertheless, the effects of the ISM are also worth considering. 
Since the majority of PNs are distributed near the Galactic plane where the mean density 
of the ISM is higher, significant interaction between PNs and the ISM is expected. Such 
phenomena can be seen from the disturbed emission on the outskirts of some nebulae (e.g., 
Xilouris et al. 1996) and from some apparently off-centered nuclei (e.g., Sh 2-216, Tweedy 
et al. 1995; KFR 1 and Me We 1-4, Rauch et al. 2000). 

The electron densities of these five PNs are all found to be higher near their equators 
than near their poles. The small pole-to-equator density ratios ((3) of some nebulae indicate 
they would show a bipolar morphology if seen edge-on. The fraction of intrinsic bipolar PNs 
therefore are likely to be higher than suggested by classification by apparent morphology. 
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Since the morphology of PNs has been proposed to relate to their physical properties, (e.g., 
bipolar nebulae have been suggested to correlate with nebulae with more massive central 
stars; Stanghellini, Corradi, & Schwarz 1993), it is therefore important that we have a good 
understanding of the intrinsic strucures of PNs. By using models such as DES, we can more 
reliably estimate the fraction of bipolar nebulae and therefore assess any correlation with 
physical properties of the PNs. 



We have derived the dust distribution in five compact PNs from their Ha and radio 
images. From this study, we are able to determine the extent of contribution by dust ex- 
tinction to the apparent departure from axi-symmetry in these PNs. By comparing these 
results with simulations generated by the ellipsoidal shell model, we can infer the intrinsic 
morphologies of the PNs. This method gives a more quantitative approach to the analysis 
of PN morphology than visual inspection of apparent morphology, and is therefore more 
capable to deducing proper statistics on the distribution of morphological types of PNs. 

We thank Kate Y.L. Su, Kevin Volk, Orla Aaquist, Steven Gibson for helpful discussions. 
Some figures were made with the software package WIP (Morgan 1995). We are grateful to 
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(Gooch 1996) 1 . This work was supported in part by grants to SK from the Natural Sciences 
and Engineering Research Council of Canada and by the National Science Council of Taiwan. 
SK acknowledges the award of a Killam Fellowship from the Canada Council for the Arts. 
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Fig. 1.— The SED of BD+30°3639, IC 5117, M 1-61, M 2-43, and M 3-45. The asterisk are 
ground-based optical and 2MASS measurements, the triangles are IRAS measurements, the 
crosses are MSX measurements, diamond are radio measurements, and the solid line is the 
IRAS LRS spectrum. The dashed line is the fitted blackbody curve for dust emission. The 
dust temperatures are BD+30°3639: T d = 190K, IC 5117: 180 K, M 1-61: 170 K, M 2-43: 
200 K, and M 3-35: 180 K. 
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Fig. 2.— Images (from top to bottom) of BD+30°3639, IC 5117, M 1-61, M 2-43, M 3- 
35. Left panel: A6 cm radio continuum (contours) on HSTHa image (color, in square-root 
intensity scale). Right panel: A6 cm radio continuum (contours) on extinction map (color). 
The color bar indicates r values. BD+30°3639: The contour levels are ax (3, 5, 10, 15) with 
the noise level a = 4.5 x 1CT 5 Jy/beam. IC 5117: the contour levels are ax (3, 5, 10, 15, 
20, 25, 30, 40, 50) with the noise level a = 8.89 x 10" 5 Jy/beam. M 1-61: the contour levels 
are ax (3, 5, 10, 15, 20, 25, 30, 35, 40, 45) with the noise level a = 8.2 x 10~ 5 Jy/beam. 
M 2-43: the contours levels are ax (3, 5, 10, 20, 30, 40, 50, 60, 70, 80) with the noise level 
a = 1.11 x 10~ 4 Jy/beam. M 3-35: the color bar indicates r values. The contour levels are 
ax (3, 5, 10, 15, 20, 25, 30, 40, 50) with the noise level a = 7.58 x 10~ 5 Jy/beam. 
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Fig. 3. — Uncertainty optical depth (contours) of (from top to bottom): BD+30°3639, 
IC 5117, M 1-61, M 2-43, M 3-35. superposed on the optical depth map in greyscale. The 
contours levels are BD+30°3639: 0.1, 0.2, 0.3, and 0.4, with 0.1 showing in white. IC 5117: 
0.02, 0.04, 0.06, 0.08, 0.1, and 0.2, with 0.02 showing in white. M 1-61: 0.03, 0.06, 0.09, 
0.12, 0.2, and 0.3, with 0.03 showing in white. M 2-43: 0.02, 0.05, 0.1, 0.2, and 0.3, with 
0.02 showing in white. M 3-35: 0.03, 0.06, 0.09, 0.12, and 0.2, with 0.03 showing in white. 
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Fig. 4. — Geometry of the dusty ellipsoidal shell model. The inclination angle i, measured 
in the ab plane, marks the direction toward the observer. 
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Fig. 5. — Model radio images assuming no radial density gradient (7 = 0). The images 
from top to bottom show the effect of inclination for % = 0°,30°,60° and 90° with different 
combinations of (3 and a. The first pair of columns has (3 = 0.9, with a = 1 and 4, 
respectively. The second and third pairs have /3 = 0.5 and 0.1, also with a — 1 and 4. All 
images are on a linear grey scale. 
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Fig. 6. — Input model dust distributions using the same parameters as in Figure 5. Each 
dust distribution simulates the dust extinction at optical wavelengths, although only the dust 
in front of the emitting gas will contribute to the optical extinction (see text for discussion). 
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Fig. 7. — Model optical images with dust extinction included. The parameters are the same 
as in Figure 5. 




Fig. 8. — Derived extinction images. The optical depths are calculated by taking the natural 
logarithm of the ratio of the model radio and optical images. The images are arranged in 
the same order as in Figure 5. 
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a) Gem Radio Continuum b) 11.7>irn Cotinuurn c) Ha Recombination Line d) Extinction Map 




Fig. 9. — Comparison of the model with the observed images for BD+30°3639. The observed 
images are in the top panels and the model images are in the bottom panels. The images 
are: a) 6-cm radio continuum emission (Bryce et al. 1997). b) Gemini 11.7/im emission (Volk 
& Kwok 2003). c) HST Ha emission (from HST archive), d) Extinction map derived from 
the radio continuum and Ha emission. All displays are on a linear scale except the HST 
observed and model images are in square-root of intensity. 
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□) 2cm Radio Continuum b) Ha Recombination Line c) Extinction Map 




Fig. 10. — Comparison of the model with the observed images for IC 5117. The observed 
images are in the top panels and the model images are in the bottom panels. The images 
are: a) A 2 cm radio continuum emission, b) HSTHa emission (Kwok, Su Sz Sahai 2003). c) 
Extinction map derived from the radio continuum and Ha emission. The radio and extinction 
images are displayed on a linear scale whereas the HST images are in square-root scale. 




Fig. 11. — Comparison of the model with the observed images for M 1-61. The images are 
arranged in the same manner as Fig. 10. 




Fig. 12. — Comparison of the model with the observed images for M 2-43. The images are 
arranged in the same manner as Fig. 10. 




Fig. 13. — Comparison of the model with the observed images for M 3-35. The images are 
arranged in the same manner as Fig. 10. 
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Table 1: Basic Information for five compact PNs 

Object PN G RA Dec Radio Diameter (") 

BD+30°3639 064.7 + 05.0 19:34:45.2 +30:30:59 6 

IC 5117 089.8 - 05.1 21:32:31.0 +44:35:48 1.5 

M 1-61 019.4 - 05.3 18:45:56.0 -14:27:42 1.8 

M 2-42 027.6 + 04.2 18:26:40.2 -02:42:58 1.5 

M 3-35 071.6 - 02.3 20:21:03.8 +32:29:22 1.5 



Table 2: The a ngular resolutions of the radio and optical images 



Object VLA HST 

Beam size (") Beam PA (°) FWHM (") 



BD+30°3639 


0.082 x 0.078 


-28.37 


0.11 


IC 5117 


0.13 x 0.12 


-85.27 


0.085 


M 1-61 


0.20 x 0.14 


15.98 


0.085 


M2 43 


0.15 x 0.13 


25.51 


0.095 


M 3-35 


0.12 x 0.12 


-81.02 


0.085 



Table 3: T he angular distances shifted in Ha images 

Object Eastward (") Northward (") Offset Distance (") 



BD+30°3639 -0.364 -0.455 0.583 

IC 5117 -0.536 -0.378 0.656 

M 1-61 -0.261 +0.261 0.369 

M 2-43 +1.548 -0.273 1.572 

M 3-35 -1.138 -0.299 1.177 



Table 4: DES mod el parameters of planetary nebulae 



Object 


a 


b 


Ab 




a 


i 


8<i> 


Se 


BD+30°3639 


18 


15 


5 


0.1 


1 


15 


0.9 


1 


IC5117 


20 


12 


5 


0.1 


1 


40 


1 


1 


Ml-61 


12 


10 


8 


0.2 


1 


70 


0.9 


1 


M2-43 


20 


12 


6 


0.3 


1 


45 


0.7 


1 


M3-35 


10 


12 


10 


0.4 


1 


80 


1 


1 



